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Abstract 

Spatial solitary waves in colloidal suspensions of spherical dielec- 
tric nanoparticles are considered. The interaction of the nanoparticles 
is modelled as a hard-sphere gas, with the Carnahan-Starling formula 
used for the gas compressibility. Semi-analytical solutions, for both one 
and two spatial dimensions, are derived using an averaged Lagrangian 
and suitable trial functions for the solitary waves. Power versus prop- 
agation constant curves and neutral stability curves are obtained for 
both cases, which illustrate that multiple solution branches occur for 
both the one and two dimensional geometries. For the one-dimensional 
case it is found that three solution branches (with a bistable regime) 
occur, while for the two-dimensional case two solution branches (with 
a single stable branch) occur in the limit of low background pack- 
ing fractions. For high background packing fractions the power ver- 
sus propagation constant curves are monotonic and the solitary waves 
stable for all parameter values. Comparisons are made between the 
semi-analytical and numerical solutions, with excellent comparison ob- 
tained. 

Keywords, colloid, solitary wave, modulation theory, stability, bifurcation. 
AMS (MOS) subject classification: 35, 78. 
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1 Introduction 

Over the last two decades the mechanical interaction between light and soft 
matter has received considerable attention, including the emergence of new 
tools in optics, such as optical tweezers and traps, see [H [12l |3]. In the col- 
loidal medium considered here, which is composed of a suspension of dielectric 
nanoparticles, the exceptionally high optical nonlinearity is due to the opti- 
cal gradient force changing the concentration or orientation of the colloidal 
particles. This leads to an intensity-dependent refractive index change and 
hence a mutual interaction between the colloidal particles and light. Pos- 
sible new applications in colloidal media include optical sensors or selective 
particle trapping and manipulation. Spatial solitary waves form in a col- 
loidal medium due to a balance between diffraction of the light beam and 
the nonlinear particle-light interaction. 

[171 1181 119) derived colloidal media governing equations by assuming 
that the colloidal suspension represents a hard-sphere gas. They used the 
Carnahan-Starling formula for the compressibility of the hard-sphere colloid. 
They then considered one and two dimensional colloidal equations and de- 
rived numerically exact propagation constant versus power curves. They 
showed that bistable behaviour occurs for some parameter values and con- 
sidered solitary wave interactions for solitary waves of the same power, from 
the same and different solution branches. They found dramatically different 
interaction behaviour for solitons from the same and different branches. In 
the 2-D case, only two solution branches can occur and the bistable behaviour 
of the 1-D colloidal solitary wave is absent. 

[20] considered a colloidal suspension of two different species of nanopar- 
ticles, one with refractive index higher than the background medium and the 
other with refractive index lower than the background medium. One species 
is described by a hard-sphere gas with the Carnahan-Starling compressibility 
formula and the other by an ideal gas. Numerical solitary wave solutions 
were found for the governing equations which show that bistability can oc- 
cur for the two-dimensional geometry, which is not possible for suspensions 
composed of only a single nanoparticle species. 

[5] modelled a colloidal suspension of nanoparticles as an ideal gas and 
considered two cases, one for which the refractive index of the nanoparticles 
was greater than the background medium, and the other for which the re- 
fractive index is lower. The governing equation in the two cases of positive 
and negative polarizability was an NLS-type equation with exponential non- 
linearity and a form of saturable exponential nonlinearity, respectively. The 
stability of the solitary waves solutions was considered numerically in one 
and two spatial dimensions. |10) considered a power series for the compress- 
ibility of the colloid particles with the coefficients of the series found by using 
the Debye-Huckel model for the particle interactions. This interaction model 
included screening effects due to the ions in the electrolyte solution. Soli- 
tary wave stability was considered, with one-dimensional waves found to be 
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stable and two-dimensional waves unstable. [16] experimentally considered 
the nonlinear optical response of a nanoparticle suspension and found that 
a higher-order NLS equation with coefficients fitted from a two term series 
expansion for the compressibility best matched the experimental data. 

The governing equation for wave propagation in colloidal media is an 
NLS-type equation for which the form of the nonlinearity in the NLS-type 
equation depends on the assumed nature of the nanoparticle interactions. 
No exact solitary wave solutions exist for these NLS-type equations, so most 
existing work on has been numerical or based on a mix of various asymp- 
totic, approximate and numerical methods. An effective technique for de- 
riving semi-analytical solutions describing the stability and evolution of soli- 
tary waves in NLS-type systems is a variational approach. This approach 
is termed modulation theory and is based on using an averaged Lagrangian 
and suitable trial functions. It has been used to study the NLS equation [14] 
and optical media, such as nematic liquid crystals [21] - Related applications 
and problems considered for nematic media include dipole motion, boundary 
induced motion and the development of bores [TTJ [3] 0] . 

In §1 modulation theory is developed for the colloidal equations stud- 
ied here. A hard-sphere colloidal model is used with the compressibility 
described by the Carnahan-Starling formula. Semi-analytical solutions for 
both one and two-dimensional solitary waves are developed. In §2 the modu- 
lation equations are solved numerically to obtain semi-analytical power versus 
propagation constant curves and neutral stability curves. These curves illus- 
trate the multiplicity of the solitary wave solution branches and the regions 
of parameter space in which unstable solutions occur. An excellent compar- 
ison between the semi-analytical and numerical solutions is obtained. In §3 
conclusions and suggestions for future work are made. 



2 Modulation equations 

Let us consider a coherent light beam (laser light) propagating through a col- 
loidal suspension of dielectric hard spheres whose diameter is much smaller 
than the wavelength of the light and whose refractive index is slightly higher 
than that of the medium in which they are suspended. With these assump- 
tions, in non-dimensional form the equations governing the nonlinear propa- 
gation of the beam through the colloidal suspension are (see [T71 ITS] ) 

du 1 

^ + ^V 2 u+(?7-7/o)w-0, M 2 ^g( V ) -g , (1) 
oz 2 

with g(rj) = ^ _ ^ + In?/, g = 0(770). 

Here u is the envelope of the electric field of the light and 77 is the packing 
fraction of the colloid particles, with 770 the background fraction. Losses due 
to Rayleigh scattering have been neglected as these losses are small in the 
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limit of the particle diameter being much smaller than the wavelength of the 
light 18 . The Carnahan-Starling compressibility approximation has been 
used. Alternative models for the compressibility alter the form of g in ([TJ. 
The Carnahan-Starling approximation is valid up to the solid-fluid transition, 
which occurs at rj = ^/2tt/9 « 0.496 in a hard-sphere fluid, see |13) . 
The colloid equations ([1]) have the Lagrangian formulation 



Here the asterisk superscript denotes the complex conjugate. The colloid 
equations (JXJ) possess solitary wave solutions. However the solitary wave solu- 
tions, in both 1-D and 2-d, have only been found numerically pTS] I19J . with no 
analytical solution known. An approximate technique which has been found 
to be useful in these cases for which there is no known analytical solution 
on which to base a perturbation analysis is the use of suitable trial functions 
in an averaged Lagrangian formulation [TT|. In the special case in which 
there is an analytical solitary wave, or periodic wave, solution this approxi- 
mate technique is the same as the modulation theory or averaged Lagrangian 
technique of Whitham J33] and other, related, perturbation techniques [14] , 
In this sense this use of trial functions in an averaged Lagrangian is an exten- 
sion of Whitham modulation theory. The approximate trial function method 
has been applied to many problems in nonlinear optics and has been found to 
give solutions in excellent agreement with numerical and experimental results 
|21[ |B] . The key is to find a good approximation to the solitary wave solution 
of a given equation. However, in many situations, the evolution of the beam, 
in particular its position and velocity, are independent, or nearly so, of the 
form of the trial function used for the solitary wave profile [23] [3] [7] . This use 
of trial functions in an averaged Lagrangian, termed modulation theory in 
analogy with the modulation theory of Whitham [23] , will be used to analyse 
solitary wave solutions of the present colloid equations. 

2.1 One spatial dimension 

Let us first consider the solitary wave solution of the 1-D form of the colloid 
equations ([TJ. The solitary wave solution of the nonlinear Schrodinger (NLS) 
equation in 1-D has a sech profile. In the small colloid concentration limit 
r\ <C 1 the colloid equations ([TJ reduce to the NLS equation. Therefore, let 
us use the trial functions 




(2) 



u = a sech 



i] = r/ + a sech' 



2 x - £ 



(3) 



w 



for the electric field and colloid fraction to analyse the evolution of an initial 
beam to a steady solitary wave solution. The parameters are all functions 
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of z. They are the amplitudes, a and a, and the widths, w and j3, of the 
electric field and colloid fraction, respectively, a is the propagation constant 
of the solitary wave, while £' = V is the velocity of the solitary wave. Lastly, 
g is the amplitude of the radiation bed on which the beam sits. The first 
term in the electric field solitary wave is a varying NLS-type solitary wave. 
The second term represents the low wavenumber (long wavelength) radiation 
that accumulates under an evolving beam which is not an exact solitary wave 
solution. The origin of this low wavenumber radiation is explained in detail 
in |14j . However, its origin can be deduced from the group velocity c g = k 
for linear waves of wavenumber k for the NLS-type equation ([I]). It can be 
seen that low wavenumber waves have low group velocity and so accumulate 
under the beam as it evolves. This flat shelf of radiation under the beam 
matches to shed radiation of non-zero wavenumber which propagates away 
from the beam and so allows it to settle to a steady state [2] . The flat shelf 
of radiation cannot extend indefinitely and so is assumed to have length £, so 
that g is non-zero in £ — £/2 < x < £ + lj2 [14]. Finally the flat shelf under 
the beam is tt/2 out of phase with it, which accounts for the i multiplying g 
in the trial function for u [14] . As previously mentioned, the flat shelf links 
with the shed diffractive radiation which propagates away from the beam. 
The effect of this shed radiation could be included [TJ]. However, for the 
present analysis the effect of this shed radiation is not needed. 

Substituting the trial functions ([3]) into the Lagrangian ([2]) and averag- 
ing it by integrating in x over the infinite domain results in the averaged 
Lagrangian 

C = -2 (2a 2 w + ig 1 ) (a' - V? + V 2 /2) - 2-nawg' + 2-ngwd (4) 

2 a 2 

+2irgaw' - - — + 4cra 2 fti(u>,^) - /9Hi(a) - 4/36i(a) 

3 w 

+4a(3 (1 + g ) , where 

f°° , ( 2 ( 
Q,i(w,B)= / sech — sech — dC, 
Jo P w 

f 00 [4-2r; -2asech 2 C 4 - 2r? Q " 
ni(ar) = 2 / 

Jo 

01 (a) = 



(l-^o-asech 2 ^ 2 (1~%) 2 



o 



?7o ln(l H sech 2 + a sech 2 £ In (rj + a sech 2 C) 

Vo 



d(, 

1 dc. 



Taking variations of this averaged Lagrangian with respect to the solitary 
wave parameters gives the modulation equations 

^- z (2a 2 w + £g 2 ) = 0, (5) 

ir^aw^igia' -V? + V 2 /2), (6) 
dz 

dg 2a 2aa ^ N /r7 ^ 
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2 --(2n 1 -wfl lw )=0, (13) 



±(2a 2 W + £g 2 )V = 0, g = V, (9) 

4a 2 n! - /3S lQ - 4/?9 lQ + 4/3(1 + 5o ) = 0, (10) 
4aa 2 fii /3 - Si - 40i + 4a (1 + g ) = 0. (11) 

Equation ((SJ) is conservation of mass, while the first of equations © is con- 
servation of momentum, in the sense of invariances of the Lagrangian ([2]). 
see [15) . However, in the present optical context ([5]) corresponds physically 
to conservation of optical power. 

A steady solitary wave solution of (JTJ will not shed radiation, so that 
g = at the fixed point of the modulation equations. Also, at the steady state 
we can set V = £ = for convenience. Note that a simple transformation 
generates non-stationary solitary waves (with V non-zero) with an unchanged 
profile. Setting g — V — S, = Q in the modulation equations gives 

1 - 3aro(fii - wQ lw ) = 0, (12) 
, 1 a 
2w^ w 

4a 2 a(ft x - - p(aZ la - Hi) - 4/3(a0i Q - ©i) = 0, (14) 

Aaa 2 il 10 - Si - 40i + 4a (1 + g ) = 0. (15) 

These four transcendental equations in the five unknowns a, a, tu, /3 and cr' 
represent a two-parameter family of solitary waves which depend on r]Q. The 
optical power is defined by 

/oo 
\u{x)\ 2 dx. (16) 
-oo 

Using the trial function (|3|) in (fT6)) gives the power of the 1-D semi-analytical 
solitary wave as 

P = a 2 sech 2 — dx = 2a 2 w. (17) 

J-oc W 

Semi-analytical power versus propagation constant curves are described by 
the solution of ([T2l) -(fT5 |) . Stable solution branches occur for P a > (see 
[IB]). Hence solitary waves of neutral stability have the property that P„ = 0. 
Adding this condition to the equations ^2l - ^5)l gives a set of five equations 
for the five unknowns, for given r/Q. Hence, the curves of neutral stability are 
lines in the a versus 770 plane. The relevant sets of transcendental equations 
for both the power versus propagation constant curves and the lines of neutral 
stability were solved using a nonlinear equation solver from the IMSL library. 

2.2 Two spatial dimensions 

The modulation equations of the previous subsection for the evolution of a 
beam in a 1-D can be extended to a 2-D beam. In this case the appropriate 
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trial functions are 

u = asech ^ e ^+iU(x- tx )+iV(,y- iy ) +i ia+iU(x- ix )+iV(y- iy ) (lg) 
W 

7] = m + a sech 2 ^, where (f> = \J {x ~ £ x ) 2 + {y - £ y ) 2 . 

Substitution of these trial functions into the Lagrangian ^) and integrating 
in x and y from — oo to oo results in the averaged Lagrangian 

C = -2 (I 2 a 2 w 2 + kg 2 ) (a' - U£' x - V£' y + U 2 /2 + V 2 /2) (19) 
—21\aw 2 g + 2I\gw 2 a + AI\awgw' — I220 2 + 2aa 2 r2 2 (u>, 0) 
-/3 2 S 2 (a) - 2f3 2 Q 2 (a) + 2/ 2 a/3 2 (1 + go) . 

In this 2-D case the shelf of low wavenumber radiation under the beam forms 
a circle of radius £, so that g is non-zero in the circle 



< ^(x-Z x ) 2 + (y-ty) 2 < L 



(20) 



Here A = ^ 2 /2. The various integrals involved in this averaged Lagrangian 
are 



Csedi C^C = In 2, 



h = / CsechC^C = 2(7, I 
Jo 

I22 = [ 
Jo 

n 2 ( w ,0) 

S 2 (a) = 







9 9 1 1 

Csech 2 Ctanh 2 C = - In 2 + -, 
3 6 

2 C 2 C 

C sech — sech — d(, 
p w 

4 - 2% - 2a sech 2 C 4 - 2r] 



(21) 



9 2 (a) 



< 



(1 - 770 - a sech 2 C) 2 (1 - ?7o) 2 

^ 2 
r/o lnfl H sech Q 

2 



+ a sech 2 £ In (770 + a sech 2 £) ] 

where C is the Catalan constant C = 0.915965594. . . 1 . The modulation 
(variational) equations for the averaged Lagrangian (fT!?)) are 



^(/ 2 a 2 W 2 + A 5 2 )=0, 

h^-aw 2 = Ag(a' - U& - V£' y + U 2 /2 + V 2 /2), 



dg 7 22 a aa 



2/i^ = 



dz w 2 



— 2 {2Vt 2 - wVL 2w ) , 



\dz dz dz 2 2 



(22) 
(23) 
(24) 
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w z 2w z 

£(l 2 a 2 w 2 + Ag 2 )U = 0, ± ( h a 2 w 2 + Kg 2 )V 



V, 



dz dz 
dz ' dz 

2a 2 n 2 - p 2 E 2a - 2(3 2 e 2a + 2I 2 (3 2 (1 + go) = 0, 
aa 2 n 2p - /3E 2 - 2/3Q 2 + 2I 2 af3 (1 + go) = 0. 



At the steady-state we have g — V — 
equations become 



(25) 
0, (26) 

(27) 

(28) 
(29) 

U = £ = and the modulation 



I 22 - a (2Q 2 - wSl 2w ) = 0, (30) 

h<T' + ^-^(4Ct 2 -wn 2w )=0, (31) 

2a 2 n 2 - (3 2 E 2a - 2(3 2 Q 2a + 2I 2 f3 2 (1 + g ) = 0, (32) 
aa 2 n 2p - PZ 2 - 2/39 2 + 2I 2 af3 (1 + g Q ) = 0. (33) 

As in the 1-D case, equations (|3U)) - (|33I) give a semi-analytical description of 
a two parameter family of 2-D colloidal solitary waves. The optical power of 
these solitary waves is given by 

y* oo 

P= r\u{r)\ 2 dr. (34) 
Jo 

The power of a semi-analytical 2-D solitary wave is then found by substituting 
the trial function (fT8|) into ([34]) . giving 

f°° r 
P = ra 2 sech 2 - dr = a 2 w 2 In 2. (35) 



3 Results and discussion 

In this section the semi-analytical solutions for colloidal solitary waves are 
compared with numerical solutions. Semi-analytical estimates for the power 
versus propagation constant and neutral stability curves are both found. 
The numerical solutions in 1-D are obtained by analytically integrating the 
steady-state version of the governing equation to obtain an energy conserva- 
tion law. The energy conservation law can then be numerically integrated to 
obtain exact solitary wave profiles on both the stable and unstable solution 
branches or the power versus propagation constant curves. In 2-D the nu- 
merical solutions were obtained using the imaginary time iterative method 
suitable for obtaining numerically exact solitary wave profiles [23] . Note that 
the imaginary time method cannot be used to obtain numerical solutions on 
an unstable solution branch. 
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Figure 1: The power versus propagation constant (P versus log a) curve for 
the 1-D colloidal solitary wave. Shown are the semi-analytical (solid lines, 
red) and numerical (dashed lines, green) solutions. The background fraction 
is 770 = 1 x 1CT 3 . 

3.1 One spatial dimension 

Figure [T] shows the power versus propagation constant (P versus logo - ) curve 
for the 1-D colloidal solitary wave. The background packing fraction is 770 = 
1 x 10~ 3 . Shown are the semi-analytical solution (fT2)) -(fT5 j) and the numerical 
solution. The figure uses the same parameters as Figure 3 of [18]. It can be 
seen that there is an excellent comparison between the semi-analytical and 
numerical solutions. The figure shows two stable branches, separated by a 
middle, unstable solution branch. The solitary waves on the low power stable 
branch have low amplitudes and large widths, whilst those on the high power 
stable branch have large amplitudes and smaller widths. The semi-analytical 
solution indicates that the middle unstable branch exists for 

-2.77 < logcr < -1.25 and 33.26 < P < 51.37 (36) 
1.31 < a < 2.29 and 4.29 x 1(T 3 < a < 0.19. 

Figure [5] shows the power versus propagation constant (P versus log a) 
curve for the 1-D colloidal solitary wave for the background packing fraction 
r/o = 1 x 10~ 2 . Shown are the semi-analytical solution (fT^]) — (fTS j) and the 
numerical solution. In this example the background packing fraction has been 
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Figure 2: The power versus propagation constant (P versus log a) curve for 
the 1-D colloidal solitary wave. Shown are the semi-analytical (solid lines, 
red) and numerical (dashed lines, green) solutions. The background fraction 
is 770 = 1 x 10~ 2 . 



increased, which eliminates the bistability. The propagation constant versus 
power curve now has a single, stable solution branch. Again the comparison 
between the semi- analytical and numerical solutions is excellent. Figures [I] 
and [2] illustrate that the semi-analytical solution is extremely accurate and 
hence is highly suitable for obtaining accurate results relating to the stability 
and other properties of 1-D colloidal solitary waves. 

Figure [3] shows the neutral stability curve in the propagation constant- 
background packing fraction plane (log cr, 770) for the 1-D colloidal solitary 
wave. Shown are the semi-analytical and numerical solutions. The semi- 
analytical and numerical neutral stability curves were found by solving the 
condition P a = using the expressions for both the semi-analytical and 
numerical power. The region under the curves represents parameter values 
corresponding to the middle, unstable branch of solitary wave solutions. The 
figure shows that as the background packing fraction increases, the region 
of parameter space in which unstable solutions occur is reduced and then 
eliminated. The parameters of the solitary wave with neutral stability at the 
turning point are 



(logCT, 770, a, a) = (-1.67, 5.69 x 1(T 3 , 1.71, 5.65 x 1(T 2 ), (37) 
(loga, 770, a, a) = (-1.69, 5.54 x 1(T 3 , 1.72, 5.79 x 10~ 2 ), (38) 
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-3.5 -3 -2.5 -2 -1.5 -1 

log o 

Figure 3: The neutral stability curve in the propagation constant-background 
packing fraction plane (log (7,770) for the 1-D colloidal solitary wave. Shown 
are the semi-analytical (solid lines, red) and numerical (dashed lines, green) 
solutions. 

where (|37|) and (|38|) are the semi-analytical and numerical solitary waves, 
respectively. It can be seen that the semi-analytical prediction of the back- 
ground packing fraction at the turning point is extremely accurate, with less 
than 3% error. Hence bistable behaviour only occurs in 1-D geometry for 
Vo < 5.54 x 10~ 3 and a single stable solution branch exists for background 
packing fractions greater than this value. This is consistent with the bistable 
curve of Figure [1] and the monotone stability seen in Figure [5] 

3.2 Two spatial dimensions 

Figure |4] shows the power versus propagation constant (P versus logcr) curve 
for the 2-D colloidal solitary wave. The background packing fraction is r/o = 
1 x 10~ 3 . Shown are the semi-analytical solution (fT2 |) -([T5 |) and the numerical 
solution. This figure is qualitatively different from the 1-D case of Figure [1] 
as the bistable behaviour is absent. There are two solution branches, one 
stable and the other unstable. The stable solution branch corresponds to 
solitary waves of large amplitude. The semi-analytical theory predicts that 
the stable branch occurs for 



a > 0.06, a > 3.81 and a > 0.46. 



(39) 
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Figure 4: The power versus propagation constant (P versus log a) curve for 
the 2-D colloidal solitary wave. Shown are the semi-analytical (solid lines, 
red) and numerical (dashed lines, green) solutions. The background fraction 
is i] = 1 x 1CT 3 . 



The comparison between the semi-analytical and numerical solutions is very 
good; there is a 17% difference in the power and a 6% difference in the 
amplitude a at a — 0.18. Note that the maximum packing fraction of the 
stable solitary wave of minimum amplitude is rj = a + t/q = 0.461, which is 
close to the value 77 « 0.496 at which solidification occurs in the hard sphere 
model. Hence this branch of stable solitary waves is unlikely to be realised 
in real colloidal suspensions. 

Figure [5] shows the power versus propagation constant (P versus log a) 
curve for the 2-D colloidal solitary wave for the background packing fraction 
r/o = 0.3. Shown are the semi-analytical solution (|30[) — (|33 |) and the numeri- 
cal solution. For this much larger value of the background packing fraction 
there is only one stable solution branch. Hence, multiple solution branches 
occur in both the 1-D and 2-D cases only if the background packing fraction 
is small enough. The comparison between the semi- analytical and numerical 
solutions is again excellent, except for some slight variation for logo - < — 3. 
For larger powers the amplitudes of the electric field and colloidal fraction 
pulses increase; beyond P ss 50 the hard-sphere colloid model predicts solid- 
ification. As in the 1-D case the comparison between the semi-analytical and 
numerical solutions is excellent and confirms the suitability of semi-analytical 
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methods for understanding colloidal wave stability. 

Figure [6] shows the neutral stability curve in the propagation constant- 
background packing fraction plane (log cr, rjo) for the 2-D colloidal solitary 
wave. Shown are the semi-analytical and numerical solutions. The semi- 
analytical solution is found by solving (f3"U |) - (|3"3"]) and the condition P a = 0. 
The numerical solutions are found using the imaginary time iterative method, 
which is used to find stable solitary wave solutions. For a given background 
packing fraction, the numerically obtained solution of minimum power is 
shown. This solution lies close to the turning point on the propagation con- 
stant versus power curve, and is a numerical estimate of the neutrally stable 
solitary wave. The region to the left of the semi-analytical curve represents 
parameter values corresponding to the unstable branch of solitary wave so- 
lutions. The curve indicates that multiple solitary wave solution branches 
only occur for 770 < 0.125 and that for larger background packing fractions a 
single stable solution branch occurs. As the semi- analytical curve of neutral 
stability is traversed from right to left the amplitude of the colloid beam 
decreases. It can be seen that there is an excellent comparison between the 
semi-analytical curve and numerical estimates for the neutrally stable soli- 
tary wave. No numerical estimates near rjg = 0.125 are presented due to 
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Figure 6: The neutral stability curve in the propagation constant-background 
packing fraction plane (log a, 770) for the 2-D colloidal solitary wave. Shown is 
the semi-analytical solution (solid line, red) and numerical solutions (circles, 
green) . 

computational difficulties in resolving marginally unstable and stable cases. 
This is related to the fact that, in this limit, the propagation constant versus 
power curve becomes very nearly flat. 

Figure [7] shows the neutral stability curve in the amplitude-maximum 
packing fraction plane (a, 77) for the 2-D colloidal solitary wave. Shown are 
the semi-analytical and numerical solutions. The semi-analytical solution 
is found by solving (|50)) - (f55]) and the condition P a = 0. The numerical 
solutions are found using the imaginary time iterative method, which is used 
to find stable solitary wave solutions. The region below the semi-analytical 
curve represents parameter values corresponding to the unstable branch of 
solitary wave solutions. This alternative view of the neutral stability curve 
shows the peak electric field intensity and colloid packing fraction of the 
neutrally stable solitary waves. As the curve is traversed from left to right 
the background packing fraction decreases from 770 = 0.125 to zero. Hence 
stable solitary waves have large maximum packing fractions when 770 is small. 
As solidification occurs for r\ > 0.496, it can be seen that physically realistic, 
but stable, solitary waves occur for a small range of amplitudes in the small 
background packing fraction limit. As for figure [51 the comparison between 
semi- analytical and numerical solutions is excellent. 
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Figure 7: The neutral stability curve in the amplitude-maximum packing 
fraction plane (a, 77) for the 2-D colloidal solitary wave. Shown is the semi- 
analytical solution (solid line, red) and numerical solutions (circles, green). 

4 Conclusions 

In this paper a semi-analytical model for colloidal solitary waves described 
by the hard-sphere gas model with the Carnahan-Starling approximation 
has been developed. Power versus propagation constant curves have been 
derived for both 1-D and 2-D geometries, with an excellent comparison with 
numerical solutions obtained. The neutral stability curves show that multiple 
solution branches only occur for low background packing fractions. The 
critical background packing fractions at which multiplicity is lost are well 
determined by the semi-analytical theory. 

Future work using the semi-analytical model will involve a number of di- 
rections. Firstly, the semi-analytical theory will be used to examine the un- 
steady evolution of colloidal beams and the development of colloidal undular 
bores. This approach has been been shown to be successful in tackling these 
problems in related optical media, such as nematic liquid crystals [3J. Sec- 
ondly, the semi-analytical tools developed here can be used to analyse other 
possible colloidal models which use different gas compressibility laws or are 
composed of multiple nanoparticle species with different optical properties 
[2Uj . This will help our understanding of the stability regimes for 2-D beams 
in experimental scenarios involving physically important colloidal materials. 
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